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Several acoustic experiments show a surprising degree of 
stability in wave fronts propagating over multi-megameter 
ranges through the ocean's sound channel despite the pres- 
ence of random-like, sound speed fluctuations. Previous works 
have pointed out the existence of chaos in simplified ray mod- 
els incorporating structure inspired by the true ocean environ- 
ment. A "predictability horizon" has been introduced beyond 
which stable wave fronts cease to exist and point-wise, de- 
tailed comparisons between even the most sophisticated mod- 
els and experiment may be limited for fundamental reasons. 
We find, by applying one of the simplified models, that for 
finite ranges, the fluctuations of the ray stabilities are very 
broad and consistent with lognormal densities. A fraction of 
the rays retain a much more stable character than the typical 
ray. This may be one of several possible mechanisms leading 
to greater than anticipated sound field stability. The lognor- 
mal ray stability density may underlie the recent, experimen- 
tally determined, lognormal density of wave field intensities 
[J. Acoust. Soc. Am. 105, 3202-3218 (1999)]. 



I. INTRODUCTION 

There is a great deal of experimental and theoretical 
interest in long-range, low-frequency acoustic pulse prop- 
agation through the deep ocean's sound channel. It has 
been investigated as a problem of wave propagation in 
random media (WPRM) and as a basis for tomog- 

raphy HQ]. Recent results from the Acoustic Engineer- 
ing Test (AET) as part of the Acoustic Thermometry of 
Ocean Climate (ATOC) project can be found in Colosi 
et al. H and Worcester et al. ||. One of the main 
challenges in analyzing and understanding long range 
acoustic propagation is in dealing with difficulties aris- 
ing from the ocean environment's tendency to generate 
multiple, weak, small-angle (forward) scattering 0. At 
sufficiently long ranges of propagation, the multiple scat- 
tering should effectively randomize an acoustic pulse so 
that it is very difficult to deduce much information. How- 
ever, several long range experiments have found a great 
deal of stability in the earlier portions of the received 
wave fronts in spite of the fluctuations inherent in the 
ocean environment In addition, it has been found 

that wave field intensity fluctuations at long range are 
consistent with a lognormal density || which would be 
reminiscent of earlier work in optics on WPRM |l^] , ex- 



cept that this earlier work was for the short range (weak 
focusing) regime. 

In the past 10-15 years, simplified models inspired 
by the ocean environment have been shown to possess 
chaotic ray limits pd|-p^|. Essentially simultaneously, 
there has been enormous progress in the understanding of 
chaotic systems |Q . Some of the most familiar emerg- 
ing concepts are simpler for bounded systems and are 
not easily applicable to open, scattering systems as we 
have here. However, there is an important tool which 
does straightforwardly generalize for our purposes, the 
stability analysis of the rays. Stability matrices can be 
constructed as a function of range for each ray. Their 
properties, such as the stability exponents, reveal the ba- 
sic character of the rays, and are at the foundation of the 
findings reported in this paper. 

There are several intriguing questions that arise from 
comparing the theoretical results to date regarding 
chaotic acoustic ray dynamics in the ocean and the high 
amount of stability observed in the data. The most gen- 
eral question concerns how an acoustic pulse - which at 
multi-megameter ranges extends to nearly 10 seconds in 
time and 2 km in depth - loses it's coherence from multi- 
ple forward scattering through interaction with internal 
waves and mesoscale energetics. Because refraction is ad- 
equate to explain the scattering physics [jl5| , the ray limit 
should suffice for understanding long-range propagation. 
Some manifestations of the underlying chaotic dynamics 
should be observed. 

It has been suggested that there exists a "predictability 
horizon" at the range of propagation defined by the scale 
over which chaotic dynamics develops |l6| [. Beyond this 
range, the wave fields should appear as random super- 
positions of many plane waves which would im- 
ply that acoustic field intensity fluctuations are Rayleigh 
distributed 19 2(J. Several problems crop up beyond 
the predictability horizon. It becomes increasingly dif- 
ficult to get numerically calculated rays to converge to 
true rays of the system. Worse, semiclassical approxi- 
mations (i.e. wave front reconstructions) from the rays 
might fail for fundamental reasons related to the break- 
down of stationary phase approximations, but one should 
recognize that more optimistic viewpoints exist on this 
issue [^T|-^3|. Whether or not this is true, it is currently 
not known to what extent tomographic inversions fail for 
any system beyond its predictability horizon where eigen- 
rays are proliferating exponentially fast with increasing 
range. In order to begin addressing these and related is- 
sues, we focus on the 'forward propagation problem' by 
performing a statistical analysis that should be much less 
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sensitive to the difficulties engendered by the predictabil- 
ity horizon. 

In fact, justifications for statistical laws derived by 
invoking stochastic or ergodic postulates are often ulti- 
mately founded on thepresence of fully developed chaos; 
see for example Ref. |2J]. Systems that once were ap- 
proached by stochastic methods have more recently be- 
gun to be regarded from the perspective of dynamical sys- 
tems. The two approaches mostly give consistent results, 
but there are important distinctions. Stochastic ray mod- 
eling is the traditional approach to the geometric limit 
of the problem of WPRM [p5| . This nondeterministic 
treatment leads one to pessimistic conclusions regarding 
the overall stability expected in an ocean acoustic pulse 
at sufficiently long range Jr|. By carefully defining the 
Lyapunov exponent, it turns out to be roughly half the 
value reported in Ref. [pUf . The scales relevant to the 
ocean are such that this factor two increase in an im- 
portant length scale might prove to be significant. Also, 
for this problem, the validity of the stochastic or ergodic 
assumptions deserves to be critically examined. It is not 
obvious that a dynamical systems perspective would lead 
to similar pessimistic conclusions as does the stochastic 
ray theory. For example, the predictability horizon con- 
cept that has grown out of the chaotic dynamics point 
of view does not necessarily lead to a sudden transition 
— regular behavior at short ranges, completely stochas- 
tic just beyond — and remnants of stability that violate 
assumptions of stochasticity could persist well into the 
horizon's initial onset. We anticipate several features of 
deterministic dynamics playing an interconnected role in 
this regard, but we focus on the importance of only one, 
intermittent-like dynamics. Intermittency is a common 
feature for nonintegrable dynamical systems |p6| , p7f . For 
the ray acoustics problem, intermittent-like behavior is 
evident through the appearance of rays which persist in 
remaining relatively insensitive to their initial conditions 
(also environment) for remarkably long ranges, as mea- 
sured on the inverse scale of the mean Lyapunov expo- 
nent. It might then be expected that the existence of 
intermittent-like dynamics might allow linear based to- 
mographic inversions based on acoustic ray models to be 
suitable to greater ranges than previously anticipated. 

One objective of this article is to illustrate the exis- 
tence of intermittent-like dynamics in the generic long- 
range ocean acoustics problem. This is a direct conse- 
quence of the wide variability in the eigenvalues of the 
stability matrix which is defined in Sect. Ill A. We demon- 
strate that the magnitude of its largest eigenvalue follows 
a lognormal distribution, and that the stability exponent 
follows a Gaussian distribution. Importantly, there is 
preliminary evidence suggesting that these distributions 
are robust, i.e. that they would be found in much more 
realistic, sophisticated ocean models Pq| . To be more 
explicit, if one knows the probability density of the sta- 
bility exponents, then one can determine the expected 
measure of intermittent-like rays that will persist out to 
the reception range. It follows that these rays will not re- 



quire extremely precise numerical interpolation schemes 
for quantities such as the gradient of the potential. 

The model upon which we rely in this paper is admit- 
tedly extremely simplified. However, it is not the model 
that is of concern, it is whether or not general features 
of simplified WPRM models carry over to the ocean it- 
self. If we are careful enough, the simplifications that 
we accept remove non-essential complications for uncov- 
ering the general physical features of interest, and no 
more. A follow-up study to this one is underway which 
uses a more realistic ocean sound speed model. It is im- 
portant for confirming the applicability of our results to 
long-range ocean acoustics experiments. 

The organization is as follows: in Sect. |n[ we introduce 
and motivate a simple model leading to a one-degree- 
of-freedom, non-autonomous Hamiltonian dynamical sys- 
tem for the rays. This is followed by a discussion of the 
analysis methods which are most critical for our study. 
They are based on the stability matrix and its well known 
properties. Sect. IV examines the fluctuation behavior of 



the stability exponents giving their densities as a func- 
tion of range. The proportion of intermittent-like rays is 
deduced and compared with the numerical results of the 
model. We finish with a discussion and conclusions. 



II. FROM WAVE EQUATION TO RAY MODEL 

We briefly outline the assumptions and approximations 
leading to the highly idealized ray model used in this 
paper. The primary physics we are concerned about in- 
volve refraction of acoustic energy due to volume inhomo- 
geneities in the ocean sound speed. We assume that inter- 
actions of the acoustic energy with both the surface and 
sub-bottom are negligible. For multi-megameter ranges 
of propagation in mid-latitude, deep ocean environments, 
a significant amount of acoustic energy is received that 
satisfies this assumption ||. As alluded to above, the 
necessary assumptions leading to the primary results are 
that: i) the linear, one-way Helmholtz wave equation is 
valid (the important point here is that backscattering is 
negligible), and ii) the spatial scales of the sound speed 
field are long compared to the acoustic wavelength so 
that ray theory is justified. A detailed derivation is read- 
ily available ■ We point out a priori that the coordi- 
nate system is three-dimensional Cartesian x = (r, y, z), 
with r the range from the source, y the transverse or 
cross-range coordinate, and z the depth from the surface. 
Thus, Earth curvature effects are neglected. 

The fundamental starting point is the linear acoustic 
wave equation |Q : 



1 d 2 ip(x;t) 
c 2 dt 2 



= V 2 V(x;0 , 



(1) 



where V'( x ; is the complex scalar wave function whose 
real part denotes the acoustic pressure. The sound speed 
field c can be taken as a function of space x only whereby 
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it has been assumed that the time scales for the propa- 
gation of the acoustic wave function are small compared 
to the time scale associated with the evolution of the 
sound speed field. For non-dispersive sources, the acous- 
tic group and phase speeds are equivalent, and one can 
linearly transform Eq. ([!]) from time to frequency, arriv- 
ing at a Hclmholtz equation with the magnitude of the 
wave vector defined as k = 2tt//c, where / is the continu- 
ous wave (CW) source frequency. Attenuation effects can 
of course be incorporated by modifying Ho be a com- 
plex quantity, but since we are interested in: 1) acous- 
tic energy that interacts negligibly with the ocean bot- 
tom, and 2) typical sources operating at frequencies with 
minimal volume attenuation (with center frequencies of 
about 100 Hz ||), ignoring attenuation effects seems rea- 
sonable. Also, one can similarly derive a reduced wave 
equation which includes variations in density, but we ig- 
nore this effect because it is known to be important pre- 
dominantly with acoustic energy that interacts with the 
ocean sub-bottom, which is not considered herein. 

The next assumption (which is quite a reasonable one) 
is that the strength of the sound speed fluctuations, 
whatever the physical process that produces them, are 
small. This allows one to neglect backscattered acous- 
tic energy, and admits the one-way Helmholtz wave 
equation, whereby one assumes a primary direction of 
propagation along the range. The so-called " standard 
parabolic approximation" is invoked next. This allows 
one to derive a linear partial differential wave equation 
of parabolic type for the complex envelope of ip. The 
principle assumption is that this envelope wave function 
evolves slowly on the scale of the acoustic wavelength. 
There many flavors of parabolic approximations that 
have varying degrees of phase errors in the complex wave 
function ip as compared to the one-way Helmholtz equa- 
tion Q, but we choose to use the standard parabolic 
approximation, which takes the form 



d<p(y, 



k dr 



T2Vi0(y,z;r) 



V(y,z;r)(f)(y,z;r) 



(2) 



where the transverse Laplacian is represented by = 
d 2 + d 2 , and the variable r is the range (propagation 
variable), but plays an exact analogous role to time in 
the Schrodinger equation of quantum mechanics. The 
parameter fco = 27r//co represents the reference wave 
number, and depends on the choice of a reference sound 
speed Co, which we take to be 1.5 km/s. The potential, 
V{y 1 z;r), is related to the sound speed fluctuations as 



V(y,z;r) 



1 



Co 



c(y,z;r) 



Sc{y,z;r) 



co 



(3) 



where the sound speed variations away from an average 
profile has been expressed as c(y, z; r) = c + Sc(y, z; r). 
The sound speed fluctuations refract the rays and lead to 



chaos in a deterministic, mathematically defined sense. 
Under the parabolic approximation, the basic problem 
maps precisely onto problems of quantum chaos |||] . The 
fields of long range acoustic propagation in the ocean 
and quantum chaos thus have the opportunity of cross- 
fertilization. 

Because the instability does not critically depend on 
having multiple degrees of freedom, we make a signif- 
icant, practical simplification in the model of ignoring 
the depth degree of freedom (z); see Ref. for a more 
detailed discussion of the model presented here. The sys- 
tem could be thought of as lying in the plane of the sound 
channel axis, but this is really just the generic problem 
of WPRM (see, for example, §). The gain in simplic- 
ity more than compensates for the loss of realism at this 
point as long as the main physical phenomena carry over 
to more realistic models. As was mentioned in the Intro- 
duction, preliminary evidence for our main results have 
been found in recent calculations incorporating a much 
more realistic model 
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The magnitude of the wave vector k is large enough 
that for the purposes of this study, we can focus on the 
ray limit. The rays can be generated by a system of 
Hamilton's equations 



dr 
dp 
dr 



8H(y,p; r) 
dp 

dH(y,p; r) 
dy 



(4) 



where y and p are the phase space variables cross-range 
(position) and horizontal slowness (momentum) respec- 
tively. The independent variable r denotes range. Corre- 
spondence with Eq. (|J) necessitates that the Hamiltonian 
is explicitly 



H = ^- + V(y;r) 



(5) 



The physical meaning of momentum is p = tan 9, where 9 
represents the angle a ray subtends in cross-range about 
the range axis. 

The state of the ocean is constantly changing, and 
its exact state is unknown. A statistical ansatz is thus 
fruitful for making assertions concerning its "typical" 
state. Assuming isotropy in the sound speed fluctua- 
tions in range and cross-range, the potential is taken to 
be a realization of a zero-mean, stationary, random func- 
tion. Thus a single correlation length scale L exists. The 
standard deviation is denoted by e = c^ 1 (8c 2 ) 1 / 2 , where 
(5c 2 ) 1 12 is the root -mean-square fluctuation of the sound 
speed about cq. Typical values for underwater acoustics 
are e = O(10 -3 ), and L = 0(100) km, but both e and L 
vary plus or minus an order of magnitude depending on 
what ocean structure is considered and the geographic 
location. For purposes of studying a fully defined, deter- 
ministic dynamical system, we complete the description 
of V by defining its correlation function to be Gaussian, 



3 



(V(y; r)V{y + 5y;r + Sr)) = e 2 exp [-{8y 2 + Sr 2 )/L 2 ] 



(6) 

We exploit this single scale throughout the rest of this 
article by transforming space variables as r — > r/L and 
y — > y/L, so the physical dimensions will always be in 
units of L. One should envision the potential as be- 
ing deterministic, even though it is selected from an en- 
semble of realizations. This implies that the potential 
is to be considered a highly complicated (albeit smooth 
and fixed) function of both y and r. To provide some 
idea of the character of this potential, contours of sound 
speed fluctuations based on a typical region of V(y; r) 
is shown contoured in Fig. |l|. The boundary conditions 
are taken as open in j/, but numerically y is treated as 
periodic, with the ray coordinate unfolded a posteriori to 
simulate the open boundary condition. A variety of ini- 
tial conditions are possible with the restriction that the 
initial momentum is always kept small enough that the 
parabolic approximation is valid all along the rays. The 
rays deriving from two such initial conditions are plotted 
in Fig. U which shows their phase space portraits (posi- 
tion, momentum) . In the absence of a varying potential, 
the solutions to the equations of motion are p(t) = po, 
q(t) = pot + qo. In this figure, rays would trace out verti- 
cal lines except in the case, po = 0, in which rays would 
show up as points, (q(t),p(t)) = (qo,0). With the po- 
tential included, the rays trace out a random-walk-likc 
motion with some drift as they move further away from 
the p — line. 
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FIG. 1. A portion of a realization of V(y; r) whose full 
domain is 20L x 320L in y and r respectively. The realization 
is constructed by the method described in Ref. 16. Contours 
are labeled in units of sound speed (m/s). The heavy solid 
contour line indicates the reference sound speed of 1500 m/s. 
The normalized root-mean-square fluctuations is e = 5x 10~ 3 . 
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FIG. 2. Two distinct ray trajectories which have travelled 
through the sound speed field characterized in Fig. |l| are 
shown out to the range of 160 L. The trajectory which starts 
at zero position and zero momentum (dashed line) is highly 
unstable, and the trajectory which starts at unit position and 
zero momentum (solid line) is stable. 

III. ANALYSIS METHODS 

The standard analysis of ray stability in the theory 
of dynamical systems begins with the stability matrix. 
From here, it is possible to calculate whether a ray is 
stable or unstable, what its Lyapunov exponent is, and 
for the unstable ray, determine the orientations of the 
associated stable and unstable manifolds that character- 
ize the exponential sensitivity to initial conditions. All 
of our results and conclusions are based on the behavior 
of the stability matrices of the rays in the model intro- 
duced in the previous section. The stability matrix is 
a strictly local analysis in range about some particular 
reference trajectory. It may be stable at one range, yet 
for a greater range be unstable. There is no restriction 
that various portions of its full history cannot have com- 
pletely different stability properties. In fact, one expects 
the portions to be almost entirely uncorrelated [ p4| . 

Often research done in chaotic dynamics uses either 
time (range) independent or periodic Hamiltonians, and 
the stability matrix is investigated about periodic orbits. 
The Hamiltonian of Eq. (S) is aperiodic, and as many oth- 
ers have done before, we slightly generalize those treat- 
ments by considering arbitrary, aperiodic rays. 

A. Stability Matrix 

The stability matrix for a ray describes the behavior of 
other rays that remain within its infinitesimal neighbor- 
hood, {5y,6p}, for all ranges. It is derived by linearizing 
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the dynamics locally; see Ref. 
the range r, one has 



321 for more details. At 
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with the stability matrix being given by the partial 
derivatives 



M = 



mil rail 
mix m22 







dp r 
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(8) 



The multi-dimensional generalizations are immediate. 
The 77i2i matrix element is well known for its appear- 
ance in the prefactor of the standard time (range) Green's 
function of the parabolic equation; it therefore gives di- 
rectly information on wave amplitudes. 
The evolution of M is governed by 



dr 



-M = KM 



(9) 



with the initial condition M (r = 0) being the identity 
matrix, and 



K = 



dydp 
d 2 H 




- 



d 2 V 
dy 2 




(10) 



The latter form is the simplification relevant for Hamil- 
tonians of the so-called mechanical type as in Eq. (||). 
Since Eq. (||) represents linear, coupled, first-order differ- 
ential equations, the elements of M can be numerically 
calculated as a function of range simultaneously with the 
calculation of its reference ray using identical numerical 
techniques, e.g. variable step, fourth-order Runge-Kutta. 



B. Stability and Lyapunov exponents 

The stability matrix has several important properties. 
It can be viewed as generating a linear, canonical trans- 
formation, and therefore its determinant is equal to unity. 
It is diagonalized by a linear, similarity transformation 



A = LML~ 



A 
A- 1 



(11) 



where the last form applies specifically to the case of a 
single degree of freedom. Here, the second eigenvalue 
must be the inverse of the first in order for det[M] = 1. 
The diagonalizing similarity transformation leaves the 
sum of the diagonal elements (trace), Tr(Af), invariant. 
It is then clear that Tr(M) is real, and three distinct 
cases may arise. The first is |Tr(M)| < 2 which is linked 
to stable motion, and it is then customary to denote 
A — exp(z6*r). The second case is |Tr(M)| = 2, and it 
is often called marginally stable because it is the bound- 
ary case between stable and unstable motion. The third 



case represents unstable motion, and is characterized by 
|Tr(M)| > 2. In Fig. §, a typically stable ray is repre- 
sented by the solid line, and a typically highly unstable 
ray is represented by the dashed line. Their distinctions 
are not immediately obvious. 

The evolution of neighboring rays about a ray that has 
|Tr(M)| < 2 satisfied from the source to the reception 
range will undergo only rotations in phase space, and sub- 
sets of phase space of finite measure where this behavior 
dominates the dynamics is precisely where intermittent- 
like rays reside. Fig. [| illustrates this characteristic be- 
havior by showing a group of stable rays winding about 
each other as they propagate. The dashed ray in the 
group is the stable ray of Fig. |^. They perform their 
"random walks", yet remain winding about each other. 
For the purposes of this paper, we make a slightly gen- 
eralized definition of intermittent-like rays as being all 
those for which |Tr(M)| remains sufficiently small over 
the range of propagation, i.e. not far from two. 
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FIG. 3. Two bundles of rays surrounding the stable and 
unstable rays of Fig. |2| (dashed lines) . For the unstable bun- 
dle, 12 rays on each side of the reference ray were chosen, each 
initially with zero momentum and uniformly sampling the ini- 
tial position over a 0.1 L window about the reference ray. The 
stable bundle used only 2 rays on each side of the reference 
ray, but used the same initial condition domain as for the un- 
stable bundle; this bundle's position is translated to 7 L. Note 
the range scale is in units of the inverse Lyapunov exponent 
v^ 1 , defined in Eq. (|l2|), and for the potential characterized 
in Fig. ^, is approximately 42.94 L. 

For unstable motion, it is customary to denote A = 
±exp(i/r) where v is positive and real. The neighbor- 
ing rays move hyperbolically relative to each other. We 
add a collection of unstable (chaotic) rays onto Fig. || to 
illustrate the distinction between neighboring groups of 
stable and unstable rays. The dashed ray is the highly 
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unstable ray from Fig. The rays were selected to span 
the same size initial neighborhood as the stable group, 
yet they fan out and become completely independent. 
For the unstable case one can introduce a definition for 
the Lyapunov exponent as 



VL 



lim 



ln(|Tr(M)|) 



(121 



Note that there is no ensemble averaging implied in the 
definition of ul. None of the theory presented thus far 
prevents it from taking on a distinct value for each ray for 
each realization of the random potential. In Sect. [V, the 



value of vl will be shown to be independent of the par- 
ticular ray, and the particular realization of the potential 
as well. It thus defines a unique length scale, v7 , which 
is used from Fig. [| onward wherever the range variable 
is involved. 

For unstable motion, |A| tends to be very large leaving 
A -1 negligible. With little inaccuracy, Tr(M) = A + 
A -1 « A , even for finite ranges. One then deduces a 
stability exponent, v from Tr(M) as 



ln|Tr(M)| 



(13) 



Thus v depends on the particular ray and varies with 
range whereas the Lyapunov exponent has no range de- 
pendence by definition. 

For any fixed range, an ensemble of v can be created 
by considering various initial conditions (by exploiting 
the isotropic and stationary properties of V(y;r)), and 
different realizations of V(y; r). The resulting statistical 
densities of |Tr(M)|, P|Tr(M)| anc ^ similarly v, p v (x), 
are the main objects of concern; the two densities are 
directly tied to each other. The cumulative probability 
distribution is given as 



F v {x) 



dx'p u (x') 



(14) 



which provides a useful tool for numerically studying the 
behavior of p v (x). It also has the utility of directly giving 
the proportion of nearly stable rays up to some argument 
set to a maximum instability criterium, v = x. 
We denote the mean and variance respectively as 



u = (v) = / dx xp v {x) , 



al = ((y - v ) ) = / dx (x- u ) p v {x) 



(15) 



where the brackets () denote ensemble averaging. For 
any real 7, ensemble averages of powers of |Tr(m)| are 
expressed as 

(|Tr(Af)| 7 ) = (exp( 7I /r)) = f dx exp(jxr)p v (x) . (16) 



Note that the case of 7 = —1/2 relates to wave amplitude 
statistics resulting from a semiclassical reconstruction of 
the wave field, and will be discussed in a future work. 



To continue the theoretical development, it is useful to 
introduce a slightly modified stability exponent, v: 



_ ln(|Tr(Af)| : 
= 27 



(17) 



Clearly v is necessarily greater than vq because of the 
important distinction of ensemble averaging before tak- 
ing the natural logarithm as opposed to the inverse order 
and the root mean square fluctuation contributions. It 
is shown in the next section that the Lyapunov expo- 
nent becomes vl = lim,.—,.^ Uq, and not lim r _ >00 v which 
surprisingly remains greater than v^. Near a parameter 
regime motivated by the ocean, we find numerically that 
analytical estimates of as being equal to v are roughly 
double their actual values. 



C. Stochastic analysis results 

An analytic estimate of v can be derived from previous 
analytic results based on a stochastic analysis which in- 
volves a strong Markovian assumption p3| , P^ , ^q| . It was 
verified in Ref. |ll| that, in the context of the present 
acoustic ray model, the stochastic analysis predictions for 
v (actually v' , see text ahead) matched to a high degree 
of precision with numerical tests. In fact, no statistically 
significant deviations were observed. Thus, although the 
stochastic system is not strictly mathematically equiva- 
lent to the deterministic dynamics, we accept the appli- 
cability of those specific results at sufficiently long ranges 
(defining this range scale is admittedly not as trivial to 
determine for the general ocean acoustics scenario as it 
is for the idealized problem). We begin with 



([Tr(Af)] 2 ) = (m 
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2(TOiim 22 ) 



(18) 



By appealing to stochastic integration techniques [|35J- 
it has been shown that in the small-e, large-r limit that 

©HP 



where 




■ exp(2z/r) , 



d 2 V(y;r-0 



(19) 
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d 2 V(y;r) 



dy 2 
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y=vo 
p=po ■ 



1/3 



= (3v^) 1/3 e 



1/3,2/3 



(20) 



(in dimensional units v' = (30r) 1/3 e 2/3 /£)- The last 
result of Eq. (^0|) is for the specific example of a Gaus- 
sian single scale potential of Eq. (||). The first result 
of Eq. ( p0| ) is more general, but requires numerical con- 
firmation for models with greater realism, and will also 



G 



depend on the ray's initial conditions for models with a 
nonuniform background sound speed field. 

By symmetry considerations of the stochastic equa- 
tions, (mfi) = (m^)- It ^ s a ^ so deduced that (7711177122) 
can, at most, grow on the same scale. Defining a corre- 
lation coefficient, 



(mu 77122) 



'22/ 



where \n\ < 1, it follows that 

([Tr(M)] 2 ) = ^(l + M )exp(27/r) 



(21) 



(22) 



Then, by using the definition of Eq. (17), one obtains 

"2 
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2r 



■In 



r(l + M) 



(23) 



Note that the second term disappears if fx equals 1/2; 
i.e. (m 11 m 2 2) = (m\ 2 )/2. We give its value numerically 
in the next section. In that case, v — v' at finite range, 
and we have an analytic estimate for v (which has not 
been derived previously to our knowledge). It is also 
worth remarking that v' is not the Lyapunov exponent 
itself (just as v is not), but rather only an upper bound. 
By analogy with the behavior of v stated at the end of 
the last subsection, v 1 will turn out numerically to be 
about double the actual v^. 



IV. FLUCTUATIONS 

The ocean is not infinite in extent, and so the dis- 
tribution of the stability exponents, v (or |Tr(Af)|) at 
a specified range r, is more directly relevant to the 
ocean acoustics problem than the Lyapunov exponent, 
Vl (or exp(^r)). In order to visualize the magnitude 
of the fluctuations we are discussing, Fig. displays the 
In |Tr(M)| for eight of the rays from Fig. gas a function 
of range out to 7.5v^ . By the right end of the figure, 
for any fixed range one ray might have a |Tr(M)| = e 13 , 
and another one might have |Tr(M)| = e°. At 7.5z^ 1 , 
there exist fluctuations in the stabilities of at least six or- 
ders of magnitude which is characteristic of broad tailed 
densities. 

To characterize the fluctuations more quantitatively, 
we consider the cumulative densities for v and |Tr(A/)|. 
An initial working hypothesis might be to check whether 
at some long, fixed range, a diagonal element of M , say 
ma, is distributed as a Gaussian random variable across 
the ensemble of V(y\ r) and derive the implied cumula- 
tive densities from there. However, there ought to be an 
identifiable mechanism for a central limit theorem (GUI) 
to be operating with respect to ma. From Eq. (0), one 
can deduce that M can be decomposed into a product of 
shorter range stability matrices. For very long r, consider 
a range Ar which is short compared to the final range r, 



yet long compared with v^ 1 . Let NAr = r where N is 
large. Then it follows that the stability matrix is given 
by the left-ordered product 



N 



M 



(24) 



1=1 



where Mi is the stability matrix for the range lAr to 
range (I — l)Ar. To a great degree of accuracy the set 
of Mi should behave independently with the only cor- 
relations being amongst the matrix elements necessary 
for maintaining unit determinant. The stability matrix 
should have the statistical properties of an ensemble of 
products of uncorrelated, random matrices. 
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1 2 3 4 5 6 7 

Range (v L ') 

FIG. 4. The stability, ln\Tr(M)\, for the reference rays 
of Fig. |^ (dashed), and 3 of their neighboring rays ini- 
tially with zero momentum and initial position shifted 
0.008, 0.016, and 0.024 L away from the reference rays. 

If there exists a limiting form for a distribution at long 
range r, one would expect the same form (with different 
parameters, i.e. mean, variance) at r/2. In other words, 
the limiting form would have to be invariant under the 
matrix multiplication process. Denoting muj as the ma- 
trix elements of Mi, for the N = 2 case, we have 



mn = W2.11TO1 11 + m 2 .i2mi,2i 



(25) 



If the mi t ij behave as independent, random Gaussian 
variables, then mn could not be Gaussian because of the 
product form. The applicability of a CLT results from 
an additive process involving random variables. Instead, 
we anticipate something closer to a lognormal density 
because the log of a product of random variables acts 
like a sum of random variables. It should be mentioned 
here that this concept has been in use in many problems 
involving statistical physics |24]. 
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To test whether |Tr(M)| is lognormally distributed, we 
calculate 100,000 rays through 5 realizations of V(y;r) 
to 7.5^7 1 (320L) (a reasonable upper bound for global 
acoustic propagation) for values of e = 2 x 10 -3 , 3 x 
10~ 3 , and 5x 10~ 3 . If |Tr(M)| is distributed lognormally, 
then v is distributed in a Gaussian manner by definition. 
In Fig. ^[ we plot the cumulative density for v. The cor- 
responding analytic Gaussian form is superposed. It is 
impossible to distinguish the numerical results from the 
Gaussian form from this plot. A similar plot for |Tr(M)| 
carries little new information, and is not pictured in this 
paper, though we have verified its excellent consistency 
with a lognormal density as well. By plotting the dif- 
ferences between the numerical and analytical curves for 
three different ranges, we see in Fig. || that the consis- 
tency with a Gaussian density is excellent, and that as 
range increases the consistency of p v (x) with a Gaussian 
density improves. Note the small scale of the deviations. 
We have verified that they are roughly of the order of 
expected sample size errors for the curve at maximum 
range. 



H? 




V (V L ) 




FIG. 6. Plot of difference between analytical and mea- 
sured cumulative densities for v at the ranges of 1.87 (light 
gray), 3. 75 (medium gray), and 7.5 (black) u£ . The measured 
cumulative densities are taken from the same simulations that 
produce Fig. ^| The values of the free parameters uo and v 
were adjusted within their simulated standard deviations to 
minimize the maximum difference for each range shown. Note 
the scale for v has been normalized by vl- 

There are several relationships implied by the lognor- 
mal density that are straightforward to test. First, if 
we denote the variance of v as a 2 , then a relationship 
between v and v$ can be derived. With 



p v {x) 



-ln(e 2 
2r X 



exp 



2r 



exp 



da; exp(2xr) 

(x-v ) 2 
2a 2 , 



FIG. 5. Plot of the cumulative density for v at the range of Inverting this last relation for a 2 vl one obtains 



7.5 vj^ 1 . The measured cumulative density is computed from 
5 realizations of sound speed fields characterized in Fig. |l| 
(e = 5 x 10~ 3 ). It incorporates 20,000 rays per realization 
whose initial conditions uniformly sample 20 L in position 
and have zero initial momentum. Superposed is the cumula- 
tive density associated with the Gaussian density for v [see 
Eq. (^6|)] using a value of uq — 0.0232 (which is our best es- 
timate for vl derived from the same simulations) and 9 = v' 
[see Eq. (po|)]. Note the scale for v has been normalized by 
vl- 



(26) 



(27) 



Both exponents D, v$ were defined (see Eq. ( |13| , |15| , |17D ) to 
be independent of r to leading order; see the upper panel 
of Fig. 7| where v, vq are plotted as a function of range. 
The stochastic approximation for v 1 also given matches 
precisely the value of v implying that fi = 1/2. From nu- 
merical simulations, it turns out that fi is 0.466, but this 
number is poorly determined due to sample size errors. 
There is no discernible r-dependence in either D or vq 
beyond the scale at which the stochastic approximation 
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begins to work for v. They maintain a rather constant 
ratio of 2.20 amongst themselves. The lognormal |Tr(M)| 
density thus implies that the standard deviation of p v (x) 
approaches zero as r~ x / 2 . Again there is excellent con- 
sistency; see the lower panel of Fig. ^ where <7„ is plotted 
versus [(D — ^/r] 1 ^ 2 . Thus, in the limit of r — ► oo, p u (x) 
goes to a ^-density; all v converge to the single value v$. 
This value would also have to be the Lyapunov exponent 
from the definition in Eq. (|l2|) , and the Lyapunov expo- 
nent would be a constant for all trajectories independent 
of the specific realization of the potential or the initial 
conditions. It appears that the approach of vq to vl is 
so rapid as to warrant replacing vq with vl in all the for- 
mulae of this section. In fact, the lower panel of Fig. 
actually incorporates our best value for vl and not vq as 
a function of range. 
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g .04 
> .02 



Stability Exponents 
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Standard Deviation 




2 4 6 8 

Range (v L ') 

FIG. 7. Upper panel: The range dependence of v derived 
from simulations (x's) with the solid line indicating the ana- 
lytical value (= 0.0510) determined by Eq. (po|). Also shown is 
the range dependence of uo derived from simulations (o's) with 
the dashed line indicating our best estimate of (vq) — 0.0232 
derived by taking the mean of the vo values at ranges beyond 
2 u~£ . Lower panel: The standard deviation of v as a function 
of range derived from simulations (o's), and the analytical es- 
timate (solid line) from Eq. ( jj^ ) using the value (yo) for vq 
(as described in the upper panel) and Eq. ( pc] ) for v. 

A consequence of the lognormal density for |Tr(M)| is 
that the density of |Tr(Af)| 7 for any real 7 must also be 
lognormally distributed. This follows from the fact that 
71/ would be Gaussian distributed with mean 71/0 and 
variance 7 2 ct 2 , and 72/ = In |Tr(M)| 7 /r. Thus 7 enters 
as a linear scale factor in the parameters that define the 
lognormal density. It is given by 



P\Tr(M)\y( x ) 



2-Kr(y — vq) \^\x 



exp 



(In (3O/7- v r) 
2r(y ~ u ) 



x>0 



Straightforward integration gives 

(|Tr(M)| 7 ) = exp ([71* + ^{v - is )/2] r) 



(28) 



(29) 



for its ensemble averaged value. Note that the 7 = 2 
case for which the stochastic theory was worked out is 
the only one independent of vq, and thus also in the 
large-r limit. Using Eq. (^9|), a variety of estimates for 
the Lyapunov exponent can be constructed. For example, 
for r large enough 



^«^ln(|Tr(M)|)-iln(|Tr(M)| 2 ) 

as given in [ 53 . Another example would be 

^«lln(|Tr(A/)| 2 )-lln(|Tr(M)| 4 ) 
r 4r 



(30) 



(31) 



etc. 

Another interesting, rather curious consequence of the 
constant ratio of v to i>q is that v does not approach 
in the r — > 00 limit even though p v {x) — > S(x — vq). Care 
must be taken to perform the non-commuting operations 
of taking the infinite range limit and ensemble averag- 
ing in the correct order. Furthermore, the variation of 
|Tr(M)| grows without bound as a function of range r, 
in spite of the fact that all the trajectories possess equal 
stability exponents in the limit r — > 00. From Eq. ([29]), 
it follows that 



|Tr(A/)| 



") 



(32) 



where the last form applies in the large-r limit, even 
though a 2 is approaching zero. 

Finally, we point out that a lognormal density has long 
tails and, as already noted, allows for many orders of 
magnitude fluctuations in |Tr(M)|. To return to the issue 
of intermittent-like rays, at any range, all rays whose cor- 
responding |Tr(M)| are less than some 0(1) constant can 
be considered as intermittent- like. Values of e or e 2 could 
be taken as criteria, for example. The equivalent crite- 
ria expressed for the maximum of v would be 0(r _1 ). In 
the present model the proportion of intermittent-like rays 
approaches zero as r — * 00, but for finite range the pro- 
portion of intermittent rays up to some maximum value 
|Tr(M)| = x is given by the cumulative density 



F, 



Tr(M)\( x ) 




ln(x) — vqt 
y/2r(D - v ) 



(33) 



where erf(z) is the error function of argument z. With 
the replacement of vq by v^, this gives a very interest- 
ing, nontrivial connection between the Lyapunov expo- 
nent {vl), and the proportion of intermittent rays as 
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a function of range. The validity of this expression is 
verified in Fig. |^. The behavior is just as predicted. The 
small deviations seen may be indicative of some slight 
non-lognormal behavior in the lower tail, or it may just 
be due to our not using best fit values of v and vq. For 
long ranges, the proportion of intermittent-like rays de- 
creases exponentially as aor 1 / 2 exp(— bor) where ao and 
60 can be deduced from the asymptotic properties of the 
error function, and this behavior is independent of the 
precise criterium used for the maximum desired |Tr(M)|. 
As can be seen in Fig. ||, 10% of the initial ray density re- 
mains stable or nearly stable out to ranges of order 5 v~^ x . 
This 10% of the initial acoustic energy is then only lin- 
early sensitive to the fluctuating sound speed field, and 
since energy remains in coherent bundles (see Fig. |^), 
they will be expected to have a longer time coherence 
over repeated experiments as the environment evolves. 
Performing repeated experiments and applying coherent 
averaging as a filter to pick up this energy, one can imag- 
ine being able to use this apportionment of the initial 
acoustic pulse for acoustic tomography. 




2 4 6 

Range (v^ 1 ) 

FIG. 8. The range dependence of the cumulative distribu- 
tion for \Tr(M) \ evaluated at e 1 and e 2 . The measured values 
(o's) are derived from simulations (see caption of Fig. [5]) and 
the theoretical curves (solid) are derived from Eq. (p3) using 
the values of v and Vq indicated in Fig. [?]. The horizontal 
dashed line at Fjiy(M)| = 0.1 indicates the range at which 
10% of the ray density remains nearly stable. 



V. DISCUSSION AND SUMMARY 

Long-range, low-frequency sound propagation in the 
ocean has been previously investigated both as a problem 
of wave propagation through a random medium, and as 
a basis for tomography. Several outstanding quandaries 
remain that our results only begin to address: i) in the 
early arriving portion of a wave front, there seems to 



be more coherence and stability than would be expected 
from an analysis based on stochastic ray techniques com- 
mon to the subject of WPRM; ii) one expects that as one 
moves from the weak focusing to strong focusing regimes 
(roughly speaking, from short range to long range), there 
should be a transition from lognormally distributed wave 
field intensities to Rayleigh distributed ones. Data anal- 
yses suggest that the lognormal densities extend well be- 
yond the weak focusing regime, and the cross-over is not 
understood theoretically; and iii) related to the first item, 
given the presence of more stability than seems consistent 
with theoretical modeling, how valid are the underlying 
assumptions of tomographic inversions performed at long 
ranges? 

Complete solutions to these problems are well beyond 
the scope of this paper. Nevertheless, we believe that 
our results form one cornerstone for their eventual res- 
olution. As the ocean acoustic problem mainly involves 
refraction, and is in a wavelength regime that should be 
extremely well-suited to semiclassical analysis, we have 
focused our attention exclusively on a simple ray model 
inspired by the ocean. Our approach is from a dynam- 
ical systems perspective as opposed to a stochastic ray 
method. It has the advantage of being a more fundamen- 
tal starting point in the sense that a system's dynamics 
may determine where a stochastic ray method is appro- 
priate, but a stochastic ray method just presupposes a 
certain randomness that may or may not actually exist 
in the system's dynamics. 

Our main concern is the ray stability properties that 
govern wave field amplitudes in semiclassical approxima- 
tions. A follow-up study is underway to address the 
phases (classical actions and geometric indices), corre- 
lations amongst ray properties, and robustness, i.e. the 
generality and applicability to the ocean of our results. 
The stability matrix is our key analysis tool because it 
contains all the necessary information about how stable 
or unstable each ray is. The distribution of stabilities 
reflects on statistical properties arising in the study of 
WPRM whereas the existence or lack thereof of stable 
rays impacts tomographic inversion. We also note that 
studying the stability matrix has the utility of providing 
additional strong checks on one's numerical integration 
techniques. Its determinant must remain unity. 

We have carefully introduced several stability expo- 
nents depending upon whether ensemble averaging is 
taken before or after the logarithm (or at all), and 
whether the range is finite or the infinite range limit is 
taken. We have related them to the absolute value of the 
trace of the stability matrix which we have found to fluc- 
tuate to a high degree of consistency with a lognormal 
density; note that this also applies to the absolute value 
of individual matrix elements of the stability matrix. We 
have given a heuristic argument for this distribution, and 
are not aware of any known analytic derivations of this 
result. 

An important consequence follows from the appear- 
ance of lognormally distributed stabilities, or equiv- 
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alently Gaussian densities in the stability ex pon ents 
(the logarithmic variables). As shown in Sect. |v| [see 
Eq. (pq)], any power of the stability matrix trace, or in- 
dividual matrix elements, is also distributed lognormally. 
Thus, each individual contribution of an cigcnray to the 
semiclassical approximation of the Green's function has 
a magnitude fluctuating as a lognormal density. Fur- 
ther study is underway to determine theoretically the 
cross-over from lognormal wave field intensity distribu- 
tions characteristic of the weak focusing limit to Rayleigh 
densities in the strong focusing limit. It is tempting 
to extrapolate our results to compute statistically rele- 
vant quantities such as the scintillation index (normalized 
variance of intensity) by using Eq. (p9|). Although one 
can immediately deduce that the normalized variance of 
intensity due to a single ray contribution grows exponen- 
tially with range with a e-folding scale of (p — vq), one 
cannot infer anything about the scintillation index in the 
region where multipathing is important since the phase 
of each contributing ray must be incorporated into the 
calculation. Also, our work assumes one is at or beyond 
the regime of strong focusing. This is confirmed in the 
upper panel of Fig. m, where vq is seen to converge at the 
range O^^ 1 ). Since the strong Markov assumption is 
valid for this problem, this range can be shown to be of 
the order where strong focusing occurs. Thus it is erro- 
neous to compute the scintillation index from our work 
in the weak focusing regime (r <C v^ 1 ) where only one 
ray contributes to the intensity distribution. 

All rays in the model possess identical Lyapunov expo- 
nents, and the finite-range, mean stability exponent, uq, 
converges rapidly to it (see upper panel of Fig. 0) . This 
follows from the finite range stability exponent acting as 
a Gaussian random variable with a standard deviation 
shrinking with range as r -1 / 2 . This is also a consequence 
of the lognormal density, and not the single scale nature 
of the sound speed fluctuations per se. This behavior 
may be rather general (as general as the lognormal be- 
havior), and it would be interesting to verify it in more 
realistic models. It is quite unlike the e 2 / 3 -scaling law 
for the stability exponents which should only apply to 
a model with a single correlation scale in range for the 
sound speed fluctuations. 

The lognormal distribution has very broad tails. One 
typically observes stability matrix traces that fluctuate 
many orders of magnitude at a given range. Long af- 
ter the appearance of highly unstable rays as a function 
of range, some stable or nearly stable rays will still be 
present. Their proportion decays essentially exponen- 
tially with range where the parameters are uniquely fixed 
by the Lyapunov exponent and the related stability ex- 
ponent v. However, they may be tomographically invert- 
ible, and relatively speaking, more important than their 
proportion would suggest. We have pointed out the dis- 
tinctiveness of their behavior relative to unstable rays 
such as the way they twist about each other, and hang 
together as they propagate. Their collective properties 
appear to be highly correlated. 



The r -1 / 2 behavior of the standard deviation of the 
stability exponent leads to a paradoxical situation in 
which all the rays possess a unique Lyapunov exponent, 
yet the exponentiated quantity, the matrix elements or 
trace of the stability matrix possess a divergent variance 
as the range approaches infinity. This illustrates dramat- 
ically the differences arising when non-commuting opera- 
tions, i.e. ensemble averaging, taking the logarithm, and 
taking the infinite range limit, are interchanged. It is 
the stability matrix elements which are relevant to semi- 
classical approximations. So the individual terms in a 
summation over eigenrays will vary infinitely in their rel- 
ative importance. 
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